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We carry out a high-precision Monte Carlo study of the shortest-path fractal dimension dmin for 
percolation in two and three dimensions, using the Leath-Alexandrowicz method which grows a 
cluster from an active seed site. A variety of quantities are sampled as a function of the chem- 
ical distance, including the number of activated sites, a measure of the radius, and the survival 
probability. By finite-size scaling, we determine dmin = 1.130 77(2) and 1.375 6(6) in two and three 
dimensions, respectively. The result in 2D rules out the recently conjectured value dmin = 217/192 
[Phys. Rev. E 82, 020102(R) (2010)]. 

PACS numbers: 05.50.+q (lattice theory and statistics), 05.70.Jk (critical point phenomena), 64. 60. ah (per- 
colation), 64.60.F- (equilibrium properties near critical points, critical exponents) 



As a standard model of disordered system [l|, Q , perco- 
lation has been intensively studied over the last 50 years 
and applied to many other fields due to its richness in 
both mathematic and physics. The nature of phase tran- 
sition of percolation has been well established. In par- 
ticular, within the two-dimensional (2D) university class, 
there are only few critical exponents left to be expressed 
exactly, among which is the shortest-path fractal dimen- 
sion d m i n , defined by 



(1) 




where r is the Euclidean distance between two sites be- 
longing to the same cluster, and I is the shortest path. 

The shortest path I between two sites in a cluster is 
the minimum number of steps on a path of occupied 
bonds/sites in the cluster, and was first studied indepen- 
dently by several groups in the early 1980's 
length I is also called the chemical distance 
latcd quantity is the spreading dimension dt 
describes the scaling of the mass TV of a critical cluster 
within a chemical distance I as M ~ £ de , and is related to 
the fractal dimension df of the cluster by di = df/d m i n . 

In percolation, the shortest-path naturally occurs dur- 
ing epidemic growth or burning algorithms. Previous 
measurements of d m in in 2D include d m in = 1.18(4) Iff, 
1.118(15) 0, 1.15(3) [11, 1.102(13) dl, 1-132(4) [l|, 
1.130(2) 0, 1.130 6(3)jl| and 1.1303(8) [3- A sum- 
mary of the early work is given in Ref. j 1 31 ] . 

In 1984, Havlin and Nossal [3 conjectured d m in = 
d f -l/v = 91/48 - 3/4 = 55/48 = 1.145833, which was 
soon shown to be too large [3EI1- In 1987 Larsson spec- 
ulated that d m in could be 17/16 or even 1, but these are 
both excluded. In 1988 Herrmann and Stanley conjec- 
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is the "red" -bond dimension and ds is the backbone 
dimension. Using Deng, Blote and Neinhuis's result 
d B = 1.6434(2) [Hi (see also [3 [11]), we find this predic- 
tion gives d min = 2 - 1.6434(2) + 0.75 = 1.1066(2), which 
is too small. In 1989, Tzschichholz, Bunde and Havlin 
[3 considered d min = 53/48 = 1.1041666..., which is 
also below measured values. 

In 1998, Porto et al. conjectured that d min is related 
to a pair-connectivity scaling exponent <?i by d m in = 
gi + /3/v where /3 = 5/36 for 2D. However, qi was later 
shown to have the exact value gi — 25/24 [3S11, which 
implies d m i n = 55/48 = 1.145833, identical to Havlin and 
Nossal's earlier conjecture flfj| - 

In 2010, one of us (Deng) and his coauthors con- 
jectured an exact expression [3 of e? m m for the 2D 
critical and tricritical random-cluster model: d m in = 
(g + 2)/(g + 18)/32<7, where g is the Coulomb-gas cou- 
pling constant, related to the random cluster fugacity q 
by relation g = (2/ir) cos -1 (<j/2 — 1). This conjecture 
is numerically correct up to the third or fourth deci- 
mal place for all values of q studied in Ref. [3- F° r 
the q — > 1 limit — i.e., standard bond percolation — the 
predicted value d mm =217/192=l. 130208 was consistent 
with the numerical results in previous works H, @, Hoi . [T2I — 
[3- I n addition, the conjectured formula exhibits other 
good properties. It reproduces the exact results for the 
critical uniform spanning tree (q — > 0) as well as for the 
tricritical q — > Potts model; at the tricritical q —> 
point, the derivative with respect to q is also correct. 

The main goal of the present work is to carry out a 
hig h-precision Monte Carlo test of the conjecture in Ref. 
jl6| in the context of 2D percolation. A numerical esti- 
mate of d m in for 3D percolation is also provided. Some 
preliminary results of this work were reported in a recent 
paper on biased directed percolation |22J. 

We simulate bond percolation on the square and the 
simple-cubic lattice by the Leath-Alexandrowicz algo- 
rithm 0, [3 j which grows a percolation cluster starting 
from a seed site. For each neighboring edge of the seed 
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site an occupied bond is placed with occupation proba- 
bility p, and the neighboring site is activated and added 
into the growing cluster. After all the neighboring edges 
of the seed site have been visited, the growing procedure 
is continued from the newly added sites. This proceeds 
until no more new sites can be added into the cluster 
(the procedure dies out) or the initially set maximum 
time step £ max is reached. 

The above procedure is also called breadth-first 
growth, and £ is equal to the shortest-path length be- 
tween the seed site and any activated sites at time step 
£. We set £ = 1 for the beginning of the growth, and 
measure the number of activated sites N(£) as a function 
of I. In addition, we record the Euclidean distance of 
each activated site i to the seed site, and define a radius 
by 



R{1) = 



if N{£) = 
if N(£) > 1 . 
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The statistical averages, N{£) = (N(£)} and K{£) = 
(R(£)), and the associated error bars are calculated. We 
also sample the survival probability V(£) that at time 
step £, the growing procedure still survives. 
At criticality, one expects scaling behavior 



N{£) ~ £ Yn , K(£) ~ £ Yr , V(£) ~ £~ Yp 



(3) 



where critical exponents Yp, Yn, and Yr are related to 
f3, v and d mia by 
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with 7 = dv — 2(3 and d equal to the spatial dimension- 
ality. In terms of exponents of epidemic processes [24j |. 
these quantities correspond to 5 = Yp, r\ = Y/v, and 
l/z-6 = Y R . 

To eliminate the unknown non-universal constants in 
front of the scaling behaviors ([3]) , we define ratio Qo {£) = 
0(2£)/0(£) for O = Af,K and V. In the £ -> oo limit, 
one has 

Y N = log 2 (QAf), Y R = log 2 (Q TC ), Y P = - \og 2 (Q v ). 

(5) 

In 2D, one has the exactly known exponents (3 = 5/36, 
v = 4/3, and 7 = 43/18 [0, In 3D, the exact 

values are unknown, and are numerically found to be 
(3/v = 0.4774(1) and v = 0.8764(7) 

Initial estimate of d m ; n for 2D. We first carried out 
simulations at the critical point p — 1/2 for bond per- 
colation on the square lattice with time step up to 
4m = 1024 and the number of samples about 2 x 10 9 . 

The asymptotic behavior of the observables Af, 1Z and 
V is expected to follow the form 
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FIG. 1: (Color online) Plot of Vl Yp - 0.045£~ 1 versus r 1 in 
2D. The Yp value is obtained via Eq. @ by setting d m i n at 
1.1302, 1.1304, 1.1306, 1.1308, 1.1310 and 1.1312, following 
the arrow. 



where higher-order corrections are neglected and the crit- 
ical exponent Y is given by Eq. (|3]). The leading finite- 
size correction exponent is known to be y\ = —0.96(6) 
— 1 [2^|. A least-squares criterion was used to fit the data 
assuming the above form. With y\ being fixed at —1, the 
fit of V gives d min = 1.1308 ± 0.0002 and b x = 0.045(5), 
and the fit of M yields d min = 1.1308 ± 0.0002. 

As an illustration, we plot V £ Yp - 0.045£ _1 in Fig. [1] 
and N£~ Yn in Fig. [51 both vs. £ _1 , where the d m i n value is 
set at a series of values in range [1.1302, 1.1312] in steps of 
0.0002, including the above estimate <i m i n = 1.1308. The 
term — 0.045£ _1 is included in Fig.[T]to remove the overall 
slope seen in the data of V; we did not do this to the 
data (Fig. [2]), and there the slope is evident. The values 
of Yp and Yn are obtained from Eq. (J4|, using the exactly 
known values of /3 and v. Because the leading corrections 
have been subtracted in Fig. [TJ it is expected that the 
curve for the correct <i m i n value should asymptotically 
become flat and reach a constant. Figure [T] shows that as 
£ increases, the curve for the conjectured value 217/192 rj 
1.1302 is bending up while the curve for 1.1312 is bending 
down. This implies that the correct d mm value should fall 
somewhere in between. A similar behavior is seen in Fig. 

where the curve for 1.1308 is approximately straight 
while those for 1.1302 and 1.1312 are bending down and 
up, respectively. 

Further simulations for 2D. Although the conjectured 
number 217/192 seems to be ruled out by the data shown 
in Figs. [T] and [U a more careful analysis is still desirable. 
The above analysis makes an assumption that the lead- 
ing correction is governed by £~ l , but the physical origin 
of this term is unclear as the leading irrelevant thermal 
scaling field has exponent yi — — 2. It is conceivable that 
more slowly convergent corrections exist but are not de- 
tected by the simulations up to £ max =1024. In particular, 
percolation can be regarded as a special case of biased- 
directed percolation with the symmetry between spatial 
and temporal directions restored [22| . In this case, multi- 
plicative and/or additive logarithmic corrections can oc- 
cur in principle such that the scaling behavior of A/", 1Z 



3 



1.68 




1.64 



0.01 0.02 0.03 0.04 

r 1 

FIG. 2: (Color online) Plot of Mf YN versus r 1 in 2D. The 
Yzv value is obtained via Eq. Q by setting d m in at 1.1302, 
1.1304, 1.1306, 1.1308, 1.1310 and 1.1312, following the arrow. 



and V is modified as 

Oil) ~ [\og(£/£ )]y™ i Y (1 + l/Ilog^/fx)]^) , (7) 

where £q and £± are constants, and y m and y c are the asso- 
ciated correction exponents. Corrections of the log log £ 
form are also possible. We note that due to cancellation 
between nominator and denominator, the multiplicative 
logarithmic correction will not explicitly appear in ra- 
tio Qo [0=N, 1Z, V), for which the scaling behavior is 
modified as 

Qo{l) = 2 Y (l + l/[\og{l/h)]v'°) , (8) 

where y' c can be equal to y c or \y m \, depending on the 
relative amplitudes of the terms associated with them. 

To investigate this, we carried out more extensive sim- 
ulations up to I max = 16384. The number of samples was 
4.5 x 10 10 for £ < 1024, 5 x 10 9 for 1024 < I < 4096, 10 9 
for 4096 <£< 8192, and 3 x 10 8 for I > 8192. 

From the Qo data, we calculate the d m i n (£) value by 
Eqs. (0| and Table Q] displays the resulting values 
of d m i n (£) from the ratios Qn, Qn and Q-p. It can be 
clearly seen that for £ < 3072, the d m i n values that derive 
from M and 1Z increase monotonically as £ increases. For 
clarity, these data are plotted in Fig. [3] The conjecture 
dmin — 217/192 would mean that the monotonically in- 
creasing curves for d mi ^ and must bend downward 
as £ become larger, and thus a very rapid drop would 
occur near origin (l/£ — > 0) in the inset of Fig. [3J which 

seems very unlikely. The % lm data are less accurate and 
not shown in Fig. [3] 

We fit the d m i n (£) data by 

dminW =dmin+M' Sl + M~ 2 , (9) 

using a least-squares criterion. The data for small I < 
^min were gradually excluded to see how the residual 
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TABLE I: Results for d m i n from Qj^, Qiz and Q-p in 2D. 
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FIG. 3: (Color online) Plot of dmin versus I in 2D, deduced 
from QaA and Q~r_. The inset shows d m j n versus r 1 . The solid 
and dashed horizontal lines correspond to d m i n = 1.13077 
and 217/192, respectively. The red (upper) and blue (lower) 
curves are obtained from the fits. 

X 2 changes with respect to £ mm . Table [TT] lists the fit- 
ting results for d^Q, d^l and d^ n . From these fits, 
we obtain d { £0 = 1.130 77(3), d^l = 1.130 77(2) and 

d lnm = 1.130 66(15), where the estimates of the error 
margins are rather conservative. We note that the co- 
efficient hi cannot be determined well in the fits for 
Anin > 32. Thus, fits with &2 = were also carried out, 
and the results agree with our above estimates of d m i n . 

We also simulated critical site percolation on an L x L 
triangular lattice with periodic boundary conditions; this 
system is known to have zero amplitude of the leading 
irrelevant scaling field with exponent yi = — 2. A row of 
lattice sites was chosen, and all the occupied sites on this 
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TABLE II: Fitting results of d m i n in 2D, for various cutoffs 
^min- "d.o.f." stands for "degrees of freedom." 
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TABLE III: Results for d m in from Qj^, Q-r. and Q-p in 3D. 



TABLE IV: Fitting results of d min in 3D. 



row were assumed to belong to the same cluster. The 
Leath-Alexandrowicz method was then used to grow the 
cluster. The chemical radius £ of the completed cluster 
was measured. From scaling I ~ £ dmin , wc determine 
^min = 1-130 7(1), also ruling out the conjectured value. 

Results for 3D. We simulated bond percolation on the 
simple-cubic lattice at the central value of the recently 
estimated critical point p = 0.248 8118(1) [33|, which is 
slightly below the previous value of p = 0.248 812 6(5) 
[341] . The simulation was carried up to f mM = 2048, with 
the number of samples 7 x 10 9 for I < 1024 and 2 x 10 9 
for I > 1024. Analogous to the procedure on the square 

lattice, we sampled N, 1Z, V and studied the ratios Q^, 

(o) 

Qiz, Qv- The values of d nli ^ deduced from these ratios 
with P/v = 0.4774(1) are listed in Table Hill The fitting 
results are shown in Table[TV]and yield d^2 = !- 375 6 ( 6 )> 
d^l = 1.375 7(6), d<2 = 1-376 5(10), and Vl = -0.7(2). 
The data of 4^)(£) and d^i) versus i' - 7 are further 



of y x . 

In conclusion, we determined the shortest-path frac- 
tal dimension d m { n for percolation in 2D and 3D to be 
1.130 77(2) and 1.375 6(6) respectively. For the 2D value, 
we use the result which follows from 1Z(£) and has the 
smallest error bars. To our knowledge, these are the 
most accurate values currently available. The conjec- 
tured value in 2D, d m ; n = 217/192 [lj|, appears to be 
ruled out with a high probability. The 3D result repre- 
sents a substantial increase in precision over the previous 
values of 1.34(1) and 1.374(4) [1]. 
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